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Classical density functional theory (DFT) of fluids is a valuable tool to analyze inhomogeneous 
fluids. However, few numerical solution algorithms for three-dimensional systems exist. Here we 
present an efficient numerical scheme for fluids of charged, hard spheres that uses O {N log A*') 
operations and O (N) memory, where TV is the number of grid points. This system-size scaling is 
significant because of the very large A'' required for three-dimensional systems. The algorithm uses 
fast Fourier transforms (FFT) to evaluate the convolutions of the DFT Euler-Lagrange equations 
and Picard (iterative substitution) iteration with line search to solve the equations. The pros and 
cons of this FFT/Picard technique are compared to those of alternative solution methods that use 
real-space integration of the convolutions instead of FFTs and Newton iteration instead of Picard. 
For the hard-sphere DFT we use Fundamental Measure Theory. For the electrostatic DFT we 
present two algorithms. One is for the "bulk-fiuid" functional of Rosenfeld [Y. Rosenfeld. J. Chem. 
Phys. 98, 8126 (1993)] that uses O(iVlogiV) operations. The other is for the "reference fluid 
density" (RFD) functional [D. Gillespie et al., J. Phys.; Condens. Matter 14, 12129 (2002)]. This 
functional is signiflcantly more accurate than the bulk-fluid functional, but the RFD algorithm 
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requires O (N^) operations. 
PACS numbers: 

I. INTRODUCTION 

Since its inception 30 years ago (reviewed by Evans [1]), classical density functional theory (DFT) of fluids has 
developed into a fast and accurate theoretical tool to understand the fundamental physics of inhomogeneous fluids. To 
determine the structure of a fluid, DFT minimizes a free energy functional il [{pk (2?)}] by solving the Euler-Lagrange 
equations 6Cl/Spk — for the inhomogeneous density profiles pk (x) of all the particle species k. This approach has 
been used to model electrolytes, colloids, and charged polymers in confining geometries and at liquid- vapor interfaces 
(reviewed by Wu [2]). Our group has applied one dimensional DFT to biological problems involving ion channel 
permeation, successfully matching and predicting experimental data [3, 4]. 

DFT is different from direct particle simulations where the trajectories of many particles are followed over long times 
to compute averaged quantities of interest (e.g., density profiles). DFT computes these ensemble-averaged quantities 
directly. However, developing an accurate DFT is difficult and not straightforward. In fact, new, more accurate DFTs 
are still being developed for such fundamental systems as hard-sphere fluids [5, 6, 7], electrolytes [8, 9], and polymers 
[10]. 

When a functional does exist, DFT calculations are, in principle, much faster than particle simulations because 
DFT requires solving only a small set of Euler-Lagrange equations. This is especially true for systems with planar, 
spherical, or cylindrical symmetry because in many cases the Euler-Lagrange equations can be integrated analytically 
over the extra dimensions. The resulting equations have only one space variable, while particle simulations are always 
performed in three dimensions. 

In systems with little or no symmetry, however, the situation is different. Many of the DFTs for important 
systems like hard spheres [5, 6, 7, 11, 12], Lennard- Jones dispersion forces [13], and electrostatic interactions [8, 9, 11] 
require computing a significant number of convolutions. This increased computational complexity quickly increases 
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computational time. Moreover, commonly-used numerical techniques scale poorly with system size, requiring O {N"^) 
operations (where N is the number of grid points). For a complex system (e.g., in biology) that requires N > 10^ for 
sufficient spatial resolution, this can, in our experience, mean the difference between 1 week of computer time for an 
O (N^) algorithm versus 1 hour for an O (NlogN) algorithm. For this reason, the vast majority of DFT calculations 
are performed in one dimension, although there are software packages for three-dimensional system. For example, 
Tramonto Software for Nanostructured Fluids in Materials and Biology has been freely available since 2007 [14]. 

For three-dimensional DFT equations, several different methods are available to iteratively solve the equations and 
to evaluate the convolution integrals. Each choice offers different trade-offs in programming difficulty, computation 
time, memory usage, and system size scalability. For example, Newton iteration requires very few iteration steps 
compared to Picard (iterative substitution) iteration, but each Newton step generally takes significantly longer than 
a Picard step. For the convolution integrals, either fast Fourier transforms (FFTs) or real-space methods can be 
used. FFTs require a regular, evenly-spaced grid and O (NlogN) operations. On the other hand, real-space methods 
can (in principle) use an unevenly-spaced grid (giving a smaller N than required by the FFTs), but require O {N^) 
operations. The Tramonto software used Newton iteration with real- space convolution evaluation. 

In this paper, we describe a FFT-based Picard iteration method. We chose this approach for several reasons. First, 
our numerical experiments showed that Picard iteration was generally faster than Newton and that in systems with 
liquid-like concentrations Newton did not always converge. Second, we found that real-space methods are impractical 
for DFT because of the specific kernels of the convolution integrals used in DFT. These convolutions integrate the 
densities pk {x) over the interiors and surfaces of spheres (described in detail in Section II) . Neither the sphere interior 
nor surface can be represented with sufficient accuracy using real-space methods; however, they can be represented 
exactly using Fourier transforms. Lastly, our solution method requires O (NlogN) operations and O (N) memory for 
hard-sphere fluids. Therefore, it scales optimally with system size. 

Currently, this optimal scalability is for uncharged hard spheres. Electrostatics is more complicated. There are 
two kinds of electrostatic DFTs in general use, both based upon a perturbation technique. In the "bulk-fluid" (BF) 
method, the electrostatic component of the free energy functional is expanded around a BF [11], while the "reference 
fluid density" (RFD) method updates the reference fluid with information from the ionic densities pk (x) [8, 9]. The BF 
method is the most commonly used (in Tramonto, for example) and we show how to implement it with the optimal 
O (NlogN) operations and O (N) memory scaling. The BF electrostatic technique can, however, be qualitatively 
incorrect [15] (see Fig. 2-4). As we describe in Section IV C, the mathematical structure of the RFD equations is 
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fundamentally different from the convolution-based DFTs of hard spheres and the bulk-fluid electrostatics method. 
In this paper we also describe an O {N"^) operations and O {N) memory implementation of the RFD electrostatics 
method. Reducing the number of operations for the RFD electrostatics method is the subject of future work. 



II. THEORY 



The DFT Euler-Lagrange equations determine the densities Pi{x) in equilibrium in the grand canonical ensemble 
which is defined by the electrochemical potential for each ion species i in the bath, /i^"*'*. The /i^"*'', in turn, are 
determined by the bath concentrations p^"*'', detailed in Appendix A. In equilibrium, the flux density for each ion 
species is identically zero, so that 



V^i = 0, 



(1) 



constraining the electrochemical potential for each ion species fii to be a constant, . 

Here the total electrochemical potential ^i{x) is a functional of the densities Pi{x), which is divided into three parts, 
an external (ext) potential, an ideal gas portion, and an excess (ex) chemical potential: 



(2) 



The ideal gas part is given by 



(3) 



where pi represents the number density of species i, k is Boltzmann's constant, and T is the Kelvin temperature 
Moreover, pf^^ is the concentration independent part of the electrochemical potential arising from an external field 
We use this to define the problem geometry, such as a hard wall. Lastly, pf^ comes from particle interactions. Thus 
in equilibrium we have 

' ^bath _ ^ext [x)~ pf (f) 



Pj (f ) = exp 



kT 



(4) 



This paper outlines an algorithm for Eq. 4 for charged, hard spheres. 

For a system of charged hard spheres, DFT decomposes the excess chemical potential into two components, the 
hard sphere (HS) and electrostatic (ES) interactions, 



/if = pfHx) + pf^{x) 

= p^^{x) + pf^{x) + z,ecb{x) 



(5) 
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where the electrostatic component is further decomposed into a mean field contribution, arising from interactions 
between uncorrelated ions, and a screening (SC) term arising from electrostatic correlations. We define Zj to be the 
valence of species i and e the elementary charge. The mean electrostatic potential </< satisfies Poisson's equation 

- eA0 = e ^ ZjPj (x) (6) 

i 

where the dielectric coefficient e is a constant throughout the entire system. The definition of the hard-sphere and 
the screening components of /i, in terms of pi constitute the heart of the density functional theory approach and are 
discussed in detail in subsequent sections. 



III. HARD SPHERE INTERACTION 

The essential DFT-specific modeling of particle interactions is contained in the definition of the chemical potentials 
and /if In order to model the interaction of hard spheres, which defines /if we use the Fundamental Measure 
Theory (FMT) [12] developed by Rosenfeld. In FMT, a suitable basis is produced which best captures the dependence 
of the potential on the densities. These basis functions, tIq,, are obtained from averages of the densities 

n^iS) = E / P^i^'Xi^' - (7) 
where the integral is taken over all space and a € {0, 1, 2, 3, V\, V2}. The weighting fvmctions ujf are given by 

where r is the spherical radial vector. Note that the VI and V2 functions are vectors, as are the associated nyi and 
ny2 functions. If constant concentrations are used in equation (7), the "fundamental geometric measures" of the hard 
spheres (surface area, volume) are recovered. 
The HS chemical potential is given by [12] 



/xf ^ (f) = kTY,j ^ K(f' )) < {x-x') d?x'. (9) 

a 

A number of different ^Hs{na) functions have been developed [5, 6, 7, 11, 12], which have different consequences, 
most notably the equation of state for a hard sphere ffuid modeled with the DFT formalism. We have used the 



anti-symmetrized version developed by Rosenfeld et al. [16], 



^HS (naj = -noln(l - ^3) + - 



247r(l-n3)' 



1 - "3 

^^ny^\\ ^^^^ 



However, other choices for '^Hsina) do not change the numerical scheme we describe below. 

It is also important to note that the 71^ integrals (7) are, up to the sign of the argument of the weight function, 
convolutions. Since the weight functions uj'^ are either even or odd, we can always convert the integral to a proper 
convolution. Therefore, they may be evaluated using the Fourier transform and the convolution theorem: 

i 

= T-^ {p. ■ cJf ) , (11) 

where J- is the Fourier transform operator and the hat denotes the Fourier image of the function. The chemical 
potential uf^^ can be calculated in exactly the same way, with pi replaced by . 

OUa 

In order to evaluate equation (11), we use the Fast Fourier Transform (FFT) for both the transformation of pi and 
the inverse transform of the product pi ■ tof. However, the ujf are distributions, and are not easily represented on 
the rectangular grid required by the FFT. Even a very fine discretization introduces unacceptably large errors and 
destroys conservation properties of the basis (e.g., conservation of total mass). Thus, if constant concentrations are 
used in Eq. 7, the geometric measures of the sphere are not recovered with straightforward real space methods. In 
three dimensions, unlike one dimension, in our numerical experiments these errors persist no matter how fine a grid 
is used. This is a severe problem for real space methods, such as those used in Tramonto. This problem might be 
resolved through a specialized quadrature, however, the authors know of no solution yet proposed. 

Rather than attempt to discretize the weight functions on a grid, we compute the Fourier transform of each weight 
function analytically, and then evaluate them on the same mesh in Fourier space as used by the FFT. The calculations 
of the analytic Fourier transforms of are given in detail in Appendix B. This strategy allows us to calculate machine 
precision convolutions with arbitrary density fields, whereas the naive discretization of the weight functions produce 
substantial errors, often in excess of the field value itself. For example, using the convolution theorem, Eq. 11, we 
recover the geometric measures for a constant density field only when using analytic Fourier Transforms of the weight 
functions. 
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IV. ELECTROSTATICS 



Mean Field 



In order to obtain 0, we solve the Poisson equation (6), for which the source is the charge density J^i-^iPi- Since 
we have access to pi from the calculation of n^, we may solve (6) in the Fourier domain, in which the Laplacian is 
diagonal. Then the mean electrostatic potential (p can be calculated by dividing by the eigenvalues of the discrete 
Fourier transform. At grid vertex j, we have 

kj) = —r '^'^'^'^^^ ^ (12) 

O, / l-cosfcj, I l-cosfc„ , l-cosfc^ \ 
\ hi hi hi J 

where h^, hy, are the grid spacings in each direction, and kx, ky, k^ are calculated as described in Appendix B. In 
order to fully specify the potential, we choose 0(0) — 0, which is equivalent to having the additional constraint 



= 0. (13) 



V 



for a domain V. 



B. Bulk Fluid Method 



The fif^ component of Eq. 5 attempts to account for electrostatic screening interactions. In the bulk fluid (BF) 
model, it is calculated as an expansion around the bath concentration. From [8] we have 



SC „ ESMth 



\x—x'\<Rij 



where Rij = Ri + Rj, Ri is the radius of ions of species i, Apj = pj — p^"*'*, c|^' (x^x') is the two-particle direct 
correlation function, ipij (x, x') is the interaction potential of two point particles of charges z^e and zje located at x 
and x', so that [17] 

c[f (f , X') + V'y (X, X') 

(15) 

~ 87T€ y 2\i\j \X-X'\ \ 2\i\j J J 

where Xk = Rk + s with s = 55^, the screening length of the bath [18, 19]. The mean spherical approximation (MSA) 
screening parameter F is derived in [19] (see also Appendix A). 

The integral in the expansion for /if*^ is a convolution, which we also evaluate in the Fourier domain. This requires 
T{Apj), calculated using the FFT, and the transform of equation (15) which is calculated analytically below. It 



(21 

should be noted that in this model of electrostatics, transformations of the c^j + ij^ij need only be calculated once, 
since they are fixed by the problem parameters. Additionally, 



T (Apj) = T (pj - Pbath) = {Pj) - T (pbath) 



(16) 



where we have already calculated ^ (pj) for the ria calculation in Eq. (11), and T {pbath) is a constant. Thus, the 
only necessary Fourier transform each iteration is the inverse transformation. 

The accuracy of the transform of Eq. 15 is key to the convergence of the nonlinear iteration for the equilibrium 
condition. In fact, we were unable to obtain convergence when evaluating these transforms numerically using the 

(2) 

FFT, and were forced to develop analytical expressions. In order to calculate each piece of c^^ + ipij, we must take 
the Fourier transform of powers of r. The generic term has the form 



in 

T 



R 



dr r"+ sin (/cr) 



47r 



(17) 



'B{R) JO 

where k is the magnitude of k. We derive a recursive definition for the integral /„ using integration by parts: 

R 



In^ f 

Jo 



For Eq.l5, we need the terms 



dr r"+ sin(A;r) 



dr r" cos(fcr) = 




cos(fcr) 







sin(A:r) 



1 R 



H^Jn n>-l 



n < -1 



n>0 



n < 



(18) 



(19) 



I-i 
h 
h 



(1 - cos{kR)) 



R 

T 

i?2 



cos( kR) + -TT sin(fci?) 



R 2 
, cos(fci?) + 2— sm( kR) - 7^ (1 - cos(fci?)) . 
k k'^ k-^ 



(20) 
(21) 
(22) 



We also need their limits as k tends to 0, 



47r o 

hm —I -I ^ 271 R^ 

k^O k 

47r ^ AttR^ 

lim ^/i = Tri?"*. 

fc-^o k 



(23) 
(24) 
(25) 



Then we have 



a(2) 



h 



\i + Aj 

-Jo 



•i\k\ \ 2AiAj 



2 /_i 



(26) 
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C. Reference Fluid Density Method 

The Reference Fluid Density (RFD) method is an alternative to the BF method to compute /if^. As shown in [15] 
and Fig. 2-4 below, it is more accurate than the BF method. The RFD electrostatic functional is detailed in [8, 9], 
and briefly summarized here. This perturbation method approximates /if^ [{Pk (v)}] with a functional Taylor series, 
truncated after the quadratic term, expanded around a reference fluid: 



[{Pkiy)}] « pf^[{pf{y)}]-kTj2 lc^^[{pf{y)};S]Ap,{x)d'x (27) 

i 

- 'y^II ^y'^^ ' ^"^^ (""^ ^'"^ 

with 

Ap, (x) = p. {x) - pf (x) (28) 

where p™^ (x) is a given (and possibly inhomogeneous) reference density profile. By defining RFD densities to be 
the bulk densities, we recover the BF perturbation method. The RFD approach makes the reference fluid densities 
functionals of the particle densities pi {x) [9]: 

pf{y)^Pk[{p^m■,y\, (29) 
where pk is the RFD functional. In [9], it is shown that the first-order direct correlation function (DCF) is given by 



where 



(^) + E / c^j\^,S-')Ap,ix')d'x' (31) 

Apk (x) = Pk (x) - Pk (x) , (32) 

c«(x) = cW[{pfc(y)};a;], (33) 

cif (f,f') = cif [{pfc(y)};f,f']. (34) 



For the RFD functional, the densities p^ (x) must be chosen so that both the first- and second-order DCFs cf'^ 
and clj can be estimated. This is possible because the densities {pk {x)} are a mathematical construct and do not 
represent a physical fluid. The particular choice of the RFD functional we use here is that of [8], which is also discussed 
in [9]: 

P^ [{Pk (x)} ; x] = 3^ / a, (x) (x) <fx' (35) 

^-KHgfj {X) J\x'~x\<B.sc(S) 
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where the {ak} are chosen so that the fluid with densities {ak {x) pk (a?)} is charge- neutral and has the same ionic 
strength as the fluid with densities {p^ (x)} at every point x. The radius of the sphere Rsc (x) over which we average 
is the local electrostatic length scale. Specific formulas for ak (x) and Rsc {x) are given in [8, 9]. In order to estimate 
the electrostatic DCFs cf'^ {x) and c'^' (a;, x') at each point, we use a bulk formulation (specifically the MSA) at each 
point X with densities pk {x), detailed in Appendix A. 

The RFD reference density p^'^^ {x) can be rewritten as the following smoothing operation 



(36) 



3 -"-SC 

where {x) = 1 — H (x) and H is the Heaviside function [20] 



0, a; < 

Hix)={ . (37) 

1, a; > 

Eq. 36 resembles a convolution, but unfortunately the screening radius Rsc'ix) is nonconstant, and thus the convo- 
lution theorem is inapplicable. We compute Rsc using (Eq. 42 in [8]) 

where Pi{x) indicates the density of species i after we have forced the mixture be locally electroneutral, but have the 
same ionic strength. 

We can express Eq. 36 in the compact notation 

p''''^{x) = J K^{x')p{x')dx' (39) 

where the kernel K^(a;') is given by 

^ "'l-'^^-g^W) . ,40) 

3 ^SC\-^) 

Since the Fourier transform is an L2 isometry, this expression is equivalent to 

p''^^(f) = J [K^(fc)] * p(fc)dfc (41) 
where we use the hat to indicate the Fourier transform and star to indicate complex conjugation. Furthermore, we 
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can calculate the Fourier transform of our kernel analytically. We have 



K^(fc) = / K^{x')e-"''-^'dx' (42) 
iE^-^\^e'^^^-^'dx' (43) 



9{\x"\-R) 
3 



4.^3 +-'dx' (44) 

47ri?3^ '''^ / d0 / dOsine / drr^e"*'^-^" (45) 









6 / dOsine I 




Jo Jo 



(46) 



where R = Rsc{x)- This integral has been evaluated above in Section IV B, so that 



r(fc) = 3e^^"--|-p^cosfci?+ pL_3infci?| . (47) 

Thus we can calculate the action of the screening operator by performing the dot product in Eq. 39 at each vertex 
of the real space grid. This algorithm has overall complexity 0{N'^), however it is accurate to machine precision. 
Alternative schemes to accelerate the operator application will be discussed in Section VII. 

In order for this formulation to be consistent, we demand that the screening radius used to construct the reference 
fluid density p^'^-^ is identical to that given by the local MSA closure. Thus, we augment our system of equations with 

Tsc [p] {x) = Tmsa [p'''^ (p)] (x) . (48) 

Here, the Ihs of Eq. 48 indicates the value of F used to determine the reference fluid density using Eq. 36, whereas 
the rhs is calculated using Eq. A6, with p'^'^^ replacing p'"'*'' as the local equilibrium value. This equation is added to 
our global system at each vertex, producing the same number of additional equations as another ion species. 

V. DISCRETIZATION AND SOLUTION IN EQUILIBRIUM 

Problem (4) is solved on a rectangular prism domain, supporting a different system size in each Cartesian direction. 
This geometry is well supported by the PETSc DA abstraction [21], which also allows for easy parallelization. The 
grid is uniform in each direction, which allows one to compute the convolutions using Fourier transform techniques. 
PETSc supports the FFTW package [22] automatically. Periodic boundary conditions are naturally enforced by the 
FFT. 

The bath potential, /i.^"*'*, and external potential, /i^^*, are calculated just once during the problem setup. The 
geometry is deflned using external potentials, /i*^^*. The excess chemical potential is dependent on the concentration. 
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as is the electrostatic potential, so these are recalculated at each residual evaluation. Moreover, the evaluation of the 
ten '^q"'^ and the at each grid point must be done at each residual evaluation since they are also dependent upon 
Pi- 

A. Nonlinear solver 

Equation (4) is a fixed point problem for each ion species i, 

Mx) = G[{pk{x')}]. (49) 

The problem is solved using a Picard iteration, since in our experiments Newton's method was both less robust, in 
that it did not always converge, and less efficient, since it took more time when it did converge. Each new iterate p^^-* 
is generated from an initial guess p^^'' using 

= G(p(°)), (50) 

where p is understood as a vector of densities over ion species. However, with higher bath densities, it is necessary to 
use a line search during the Picard update rather than just successive substitution. Thus, our new guess p* is given 

by 

p* = (l-a)p(i)+ap(°), (51) 

where a is the line search parameter. We determine a by sampling the function G at several densities, fitting the 
residual values, \\p — G{p)\\, to a polynomial in a, and choosing amin corresponding to the minimum residual value. 
We currently have a quadratic line search, suggested to us by Roland Roth, which fits the squared L2 norms of the 
residuals from Eq. 50 as this seemed to better match curves in the search parameter we sampled for testing. 

In addition, because n3(x) is the local packing fraction, it should never exceed unity. We bound it by 0.9 which 
allows us to bound the maximum allowable search parameter a since 713 is a linear function, 

||n3((l - a)p(°)) + n3(ap(i))|U < 11(1 - a)^3(p^°))l|oo + ||an3(p«)|U 

= (l-a)||n3(p(°))|U + a||n3(p('))||oo 
< 0.9. 

Here, HafHoo is the L°° norm, which picks the maximum value of x in the finite dimensional case. This bound was 
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also suggest by Roland Roth. Finally, we have 

^ o.9-IK(p<°^)IU 

||n3(pW)||oo-||n3(p(0))||oo- ^ ' 

We have also experimented with Newton's method, forming the action of the Jacobian operator using finite differ- 
ences. The linear systems are solved with GMRES. Both fixed linear system tolerances and those chosen according to 
the Eisenstat- Walker scheme were used. However, the Newton method was not competitive with Picard due to linear 
convergence through most Newton steps and the large cost of computing the Jacobian action. 

B. Numerical Stability 

With a coarse grid, there is a potential for serious roundoff error when calculating both the average over an ion 
surface, n2, and the directional average, nv2- From the definition (7) we have: 

n2{x) = XI / P^{x')^f{^-x')d^x' (53) 

= J2 f P^i^W'^'\ - R^)d^^' (54) 

= / p^{x + ^^)dn. (55) 
Here r = i?j(sin0cos(/),sin(?sin0,cos6') and S{Ri) is the surface of a sphere of radius Ri. Likewise, 

^(^) ^ E / P^i^>Y^i^~^')d'x' (56) 

= Y.J P^i^'f-HW~£'\-R^) (57) 

i 

= 2_] / pi(x + r)(sin0cos0,sin0sin(/),cos6')(iri. (58) 



nv2[ 



Appendix C shows that 

^ < 1. (59) 
\nv2\ 

However, discretization errors in the computation of the last term of Eqn. (10), or its derivative can combine 
to produce large values of this ratio, stalling the nonlinear solve and leading to unphysical artifacts. These artifacts 
produce large density oscillations at sharp corners along the geometric boundary. These oscillations eventually cause 
divergence of the nonlinear iteration and prevent accurate solution of the equations. We alleviate this problem by 
enforcing the bound explicitly. 
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VI. VERIFICATION 



At several points in the calculation, we perform consistency checks of the results. Moreover, we compare our results 
to known thermodynamic solutions, in the limit of very fine meshes. We first check that we recover the fundamental 
measures, which are easily computed analytically, when we compute na{p) with a constant unit density. We also 
check that in the bath, 713 is equal to the combined volume fraction of the ion bath concentrations. 

Moreover, we can verify that in symmetric situations, such as near a hard wall, the solution must be homogeneous 
over each plane parallel to the wall. As described earlier, these consistency checks are satisfied when analytic Fourier 
transforms of the weight functions loi are used, but not using the FFT or real space methods. We can solve an 
effectively one-dimensional problem with a wall at z = and periodic in each dimension in order to compare with 
thermodynamic results. With purely hard sphere interactions, we have another consistency check, namely a relation 
between the pressure P in the bath and the density of each species at the wall [15, 23], 



where we evaluate the density at the point of closest approach to the wall. The FMT DFT of hard spheres is known 
to satisfy Eq. 60, but this relation holds only approximately for electrostatic functionals described here [15]. 



A very sensitive test for calculations of ionic solutions are thermodynamic sum rules, such as Eq. 60. We use this 
as the figure of merit to assess the accuracy of our hard sphere calculations. 

A notable advantage of the DFT formulation over particle simulations, such as a Monte Carlo for hard spheres, is 
that both very low and very high densities can be handled efficiently with no algorithmic changes. Low densities are 
difficult for canonical ensembles, such as canonical MC or Molecular Dynamics (MD), because very large systems are 
required for accurate statistics. A grand canonical formulation of MC can mitigate the problems for low densities, 
however, high densities still result in jamming and high rejection rates, requiring very long run times. This can 
sometimes be repaired using very specialized techniques [24, 25], however currently these cannot be applied to general 
systems of the type we present below. 

In order to demonstrate the performance of our algorithm across a range of densities, we simulate a hard sphere 
liquid against a hard wall. The particles have radius O.lnm. In Fig. 1, we show both the simulation time and accuracy 
over volume fractions ranging from to 0.4. Our results are quite accurate, and even at liquid densities the 




(60) 



A. Hard sphere fluids 
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Sum Rule Verification against Hard Wall 




packing fraction r] 



FIG. 1: Both accuracy and simulation time are shown for a hard sphere liquid of O.lnm particles. The domain is divided into 
cubes which are 0.05nm x 0.05nm x .00625nm. 

calculations done on a laptop take less than 1.5 hours. Note that, although this is an effectively one dimensional 
problem to facilitate verification, the computation was performed in a full three dimensional geometry. 

B. Ionic fluids 

Calculation of ionic densities near a hard wall also provides a sensitive test for the accuracy of the DFT method. 
In [15], it is demonstrated that the BF version of DFT (Eq. 14-15) provides qualitatively incorrect densities, when 
compared with the RFD functional (Eq. 28-35) and high resolution Monte Carlo simulation. We have successfully 
reproduced the one-dimensional DFT and Monte Carlo results with the 3D code, attesting to the accuracy of our 
approach. Below, we detail a representative simulation. 

For our trial calculation, we examine a salt solution of univalent ions. The cation has radius O.lnm, the anion 
0.2125nm. Each species has a IM bath concentration. The simulation cell, 2nmx2nmx6nm, is periodic in each 
direction. A hard, uncharged wall is placed a z = 0. We discretize the density on a 21x21x161 grid. The results 
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FIG. 2: Comparing ID and 3D DFT cation concentrations to MC simulations. The wall is uncharged, the cation concentration 
is IM, and the ions are univalent. The ID RFD DFT is shown with the blue squares, 3D RFD DFT with blue circles, and MC 
with green squares. 

are insensitive to the resolution in the transverse (x — y) directions, but very sensitive in the normal (z) direction. 
We verify the homogeneity of the solution across x — y planes to machine precision. In Fig. 2 and Fig. 3, we show 
the excellent match between ID and 3D DFT results, with MC results shown for comparison. The mean electrostatic 
potential is shown in Fig. 4, also with good agreement. 

The BF calculations are currently much more efficient than the RFD calculations, needing only 1.5 minutes com- 
pared to more than a day to run, since BF scales as 0{NlogN), whereas the RFD method scales as 0{N^) and 
additional iterates are needed to obtain a converged reference density. However, the extra investment of time for RFD 
computations is necessary because the BF solution is qualitatively incorrect compared to Monte Carlo simulations. 
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FIG. 3: Comparing ID and 3D DFT anion concentrations to MC simulations. The wall is uncharged, the cation concentration 
is IM, and the ions are univalent. The ID RFD DFT is shown with the blue squares, 3D RFD DFT with blue circles, and MC 
with green squares. 

VII. CONCLUSIONS AND FUTURE WORK 

We have presented a full numerical strategy for solving the three dimensional equilibrium DFT system. The hard 
sphere calculation accurately reproduces thermodynamic sum rules, and agrees with prior MC simulation. Moreover, 
using the improved RFD electrostatic formulation due to Gillespie et al. [8], we can accurately reproduce electrostatic 
behavior near a hard wall for species of differing radii. Thus, the DFT can now become a powerful tool for full three 
dimensional chemical simulation, accurately capturing both the energetic and entropic contributions to the solution. 

There are also several avenues for improvement of the RFD algorithm and extension of the capabilities of the current 
code. The dominant cost of this algorithm is the calculation of the reference density used to describe electrostatic 
screening. The current algorithm is very accurate, but requires 0{N'^) work. Since the Fourier kernel is smooth 
and has rapid decay, it should be possible to construct a multiresolution analysis of it, resulting in a fast method 
for application. Moreover, the many FFTs performed at each Picard step could be replaced by Unequally-Spaced 
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Electrostatic Potential 
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FIG. 4: Comparing ID and 3D DFT mean electrostatic potential to MC simulations. The wall is uncharged, the cation 
concentration is IM, and the ions are univalent. The ID RFD DFT is shown with the blue squares, 3D RFD DFT with blue 
circles, and MC with green squares. 

FFTs or wavelet decompositions, which would allow adaptive refinement and increase the size of problems we can 
efficiently compute. The FFT and Fast Wavelet Transform (FWT) lend themselves readily to a scalable parallel 
implementations. In fact, it should also be possible to offload these transforms onto a multicore co-processor, such 
as the Tesla 1060C GPU [26]. This will make large scale simulations of charged hard spheres accessible to working 
scientists even on a laptop or desktop computer. These algorithmic improvements are the focus of current research. 
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APPENDIX A: CALCULATION OF THE BATH CHEMICAL POTENTIAL 

Here we describe the formulas for the electrochemical potential in a homogeneous fluid. When the DFT for hard 
spheres uses a Percus-Yevick equation of state [27], and the electrostatics is described using MSA [18, 19]. We follow 
the treatment in [28]. The bath chemical potential /i^"*'* has two components, hard sphere and electrostatic, 

,hath _ ,HS.bath . ES,bath /' A 1 ^ 

l^i — H'i "T Mi {-^^J 

which are calculated thermodynamically [29], 

^HSM.. ^ ,T (- In A + 'J^-l±^ + + (A2) 



A 2A2 6kT 

based upon the auxiliary variables, where (7^ is the ion diameter of species i 



^» = IY.pT"'''? ne{0,...,3} (A3) 

3 

A = 1-^3 (A4) 



and the pressure due to hard sphere interaction in the bath. 



The calculation of ^^•S'.''"*'' given above is straightforward, but ^^•5'''°*''^ on the other hand, is dependent on an 
implicitly defined parameter F, the MSA inverse screening length, 

(AO) 

where 77 represents the effects of nonuniform ionic diameters 

77-- — V^^^^ iA7) 
3 ■' 



and is determined by 

„bath _3 



n 

2I\ ■ 

3 

This implicit relationship is a quartic equation in F, which we solve using Newton's method. We may then calculate 
the bath potential 



ES.bath 



l + Ta, ' V 1 + Fa, 



(A9) 
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APPENDIX B: EVALUATION OF THE FOURIER TRANSFORM OF THE WEIGHTING FUNCTIONS 

We must be careful to evaluate our analytic transforms at the same k values, in the same order, as those computed 
using the particular implementation of FFT we use. Given a D dimensional grid, the vector {kd} which corresponds 
to the vertex {j^} of our Cartesian grid is given by 

kd={ (Bl) 

-2^(jVrf-jrf) . JVi 
Naha ■'d ^ 2 

where is the number of grid points in dimension d £ {x, 2;}, and hd is the grid spacing ■ 
We begin with the calculation of cj^, 

\ d(f> dOsinO dr r^5{\r\- R{)e-'''-^ (B2) 

27r /"TT ^ 

d0 sin e ^. (B3) 

We now choose a rotated coordinate system (the prime system) in which k points purely in the z' direction, in order 
to take advantage of the rotational symmetry of the problem. In the new coordinate system, 

^1 = d<j)' de' sine'Rle~'^^^'^''°'''^' (B4) 
Jo Jo 



2ttR} [ de' sine' {cos {Rik'^ cos e') ~ ism {Rik'^ cose')) (B5) 
"'0 

AirRi sin {Rik'^) 



k'^ 

which, in the original coordinate system, is 



From Eq. (8), we also have 



(B6) 



AnRi sin ( RM] 

= ^- (B7) 



^ sin (^Ri\k\j ^ sin (^R^\k\^ 

'^i = — ;rTT^ — '^i = 7T^, ■ (^8) 



Recognizing that the theta function can be obtained as the integral of a delta function, we have 

ujf^ I drCol 1^.=, (B9) 

<^R^ 







drr sin [r\k\] (BIO) 



47r 

|fc| Jo 

p- (sin(i?,|fc|) -i?,|fc|cos(i?,|fc|)) (Bll) 
\k\ 
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Following a similar procedure as in the i2)f calculation, but keeping track of the vector nature of and w 



V2 



\V2 



Jo Jo 



2 ^—iRik'^ cos 9x 



= -2mR, / d9' sine' cose' sin {Rik'^ cos e')k 

Jo 

= (sin (Ri\k\j - Ri\k\cos (-Ri|fc|)) k. 



(B12) 
(B13) 
(B14) 



The preceding expressions for ujf may be evaluated at \k\ = 0, but care must be talcen when calculating the limit. 



sm 



lim d)^ = lim 
|fc|->o |fe|->o Ri\k\ 

lim ibl = Ri 
\k\^o 



(i?,|fc|) _ ^ 



lim ibf 

IfeHo 



lim ujf 



AnRf 



lim ^— 
|feHo|A;|3 



Ri\k\ - 



R^\k\ 



\ 



^ {R^\k\y 



lim 

IfeHo 



3 * 




Um cjf 2 = 

IfeHo 

It should be noted these are the limits one would expect, since in the |fc| = case we are simply integrating either a 
spherical delta or step function over all space, thereby recovering surfeice area and volume expressions for a sphere. 

APPENDIX C: DIRECTIONAL AVERAGE BOUND 



We can bound the directional average of the density over a sphere in terms of the unweighted average, and thus we 
can bound the ratio 



(CI) 
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in the calculation of $(n) from Eq. 10. We let ^{6, 0) be the unit vector at the surface of the sphere in the {6, t 
direction. Using Fubini's Theorem and the Cauchy-Schwarz inequality, we have 



so that 



(x) = V / iy{e,(j))p,{x + r)dn- [ v{e\(t)')pj{x + r')dn' (C2) 

= / iy{e,(f,) ■ i^{0',(l)')pi{x + r)pj{x + r')dndn' (C3) 
^ 12 [ W{0,(l))-iy{e',(j>')\\p.,{x + r)pj{x + r')\dndn' (C4) 

Js-^xS^ 

^ [ W{0,(|))Me\(|)')\p^{x + r)pJ{x + r')dndn' (C5) 

Js'^xS^ 

< V / p^{x + r)pj{x + r')dndn' (C6) 

Js'^xS^ 

= / Pz{x + r)dfl [ pj{x + r')dn' (C7) 

''J 

= nlix) (C8) 



J ^ < 1. C9 
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